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ABSTRACT 


The classical approach for the analysis of linkages has been to 
assume rigid elements during operation. As the operational speeds 
of mechanisms increase, this assumption creates larger errors in the 
analysis. Thus a better model of a machine that considers elastic 
deformations will aid in more efficient design as well as provide 
improved accuracy and performance. 

This thesis presents a simulation for the dynamic response of 
an elastic link model with three clearances. The simulation has 
been developed for an automobile, cam actuated valve train. The 
computer program is coded in FORTRAN IV, for an IBM-360 computer 
and is included as a portion of this work. The capability to visually 
observe the dynamic action of the model is included by a graphic 
display routine. This routine is implemented in SDS FORTRAN IV 
for a SDS-9300 computer Pepe caeed with an ADAGH Graphic Display 
Terminal, Model 10. 

Sample problems for various cam speeds of interest are included 


utilizing graphs and photos from the graphic display. 
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I. INTRODUCTION 


The classical approach to the dynamic analysis of linkages has 
been to assume that all menbers of the machine move without deformation 
during operation, i.e. that they move as rigid bodies. Based on this 
assumption, material choices and sizing to withstand stresses for the 
applied loads, the linkage deformations and the response motion have 
been severly idealized. Application of factors of safety have been 
necessary to account for this idealization. At lower operational 
speeds design criteria are met. However as the speeds increase the 
resulting large inertial forces and elastic deflections are more 
severe and thus the possibility of failure more acute. The design 
becomes invalid when the critical speed is reached where the machine 
either fails or where it no longer performs its intended task. Thus 
expensive and time consuming models must be built and tested, for 
elastic deformations over the range of operational speeds before the 
kinematic analysis can be fully accepted. The goal then is to obtain 
amore accurate dynamic model of a machine in order to better predict 
and understand its performance. 

With the advent of the digital computer, more comolex math= 
ematical models of machine are being developed. lEngineers are able 
to accomplish more refined designs with a higher degree of operational 
confidence, less of a factor of safety, better material choice, and 
fewer physical model experiments. 

This work presents a simulation for the dynamic response of a 


planar elastic link mechanism, free to move in its rigid body motion 
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as well as its elastic motion. The structural elements deform 
elastically due to forces applied by neighboring elements, and due 
to inertial forces. This simulation has been developed for an auto- 
mobile, cam actuated valve train, with three clearances. The clearances 
are considered to exist at the (1) cam-pushrod, (2) valve-valve seat, 
and (3) rocker arm-pushrod. To simplify the required derivations, 
the model is assumed to initially consist of straight members with 
constant circular cross-sectional bars. The assumption of constant 
circular cross-sectional bars is selected for ease in computing beam 
section properties. The use of straight links greatly simplifies the 
equations of elastic deformation. 

An existing computer program written by Anderson and Winfrey 
(Ref. ELE was dependent on separately coded sections for each physical 
configurations of the two clearance model of a cam actuated valve 
train. As the number of system clearances are increased, the number 
of system configurations increase exponentially and separately coded 
program sections become infeasible. 

The extent of this author's objectives and work on the existing 
program is to: 

(1) Restructure and reprogram the simulation to provide for a 
generalized method for establishing the physical model configurations. 
(2) Modify this program to allow for more than two clearances, 

i.e. add clearances internal to the systen. 
(3) Establish a method for the dynamic graphical display of the 
model showing the resulting elastic deformations of the model members 


at high operational speeds. 
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Objectives (1) and (2) have been accomplished using the IBM-360 
computer for the simulation. Objective (3) has been accomplished 
with the aide of the ADAGE Graphics Terminal model 10 (AGT~10) 
interfaced with a Scientific Data System computer model 9300 (Photo 
We). Data generated by the IBM computer is loaded into the SDS 
computer by means of punched cards or magnetic tape. This data is 


then manipulated into a form acceptable for display by the AGT-10. 
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In the following chapters, the physical model, and the concepts 
and formulation for the corresponding mathematical model are described. 
The dynamic analysis and the elastic deformation computer programs 
are discussed and listed at the end of this work. Finally, several 
numerical examples are given with the model responses shown as AGT 


photographs and CALCOMP graphs. 
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In the following chapters, the physical model, and the concepts 
and formulation for the corresponding mathematical model are described. 
The dynamic analysis and the elastic deformation computer programs 
are discussed and listed at the end of this work. Finally, several 
numerical examples are given with the model responses shown as AGT 


photographs and CALCOMP graphs. 
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II. PHYSICAL MODEL 


in this chapter the physical system, Figure Cols 1S. idealized to 
enable formulation of a mathematical model whose analysis is tractable. 
The description covers the interaction between the system parameters 
and the physical contraints that limit the model's motion. The 
physical elements that comprise the system are listed and their 
functions explained. Additionally their physical properties, size 


and material content, are noted. 


RIGHT ROCKER ARM 
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Figure 2.1 CAM ACTUATED VALVE TRAIN 
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A. GENERAL DESCRIPTION 
In order to relate the formulation of this mathematical model 
to the real world, an automobile car engine valve train assembly is 
considered. However, any similar physical model might have been 
chosen. All members of the model are assumed to be manufactured from 
steel with an elastic modulus of 30 x 10° pounds per square inch and 
a specific weight of .283 pounds per cubic inch. Additionally long- 
itudinal members are considered as rods with circular cross sections. 
The diameter and length of each bar are specified. 
1. Elements of the Model 
The elements of the valve assembly pictured in Figure De ieanes 
a. Follower 
This reciprocating flat faced follower is rigidly attached 
to the pushrod. The pushrod is considered as an elastic member. 
These members have harmonic displacement. The diameter is .3 inches; 
length 9 inches. 
b. Left Rocker Arm 
The left rocker arm is pinned to and rotates about the 
fixed pivot point at its right end. This element provides for an 
"tanternal" clearance between itself and the pushrod at its left end. 
The diameter is .6 inches; length 1 inch. 
c. Right Rocker Arm 
The right rocker arm is also pinned to and rotates about 
the fixed pivot point. At its right end it is pinned to the valve 


stem and valve spring. The diameter is .6 inches; length 1.5 inches. 
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a Valve Spring 
The free end of the valve spring is connected to the 
rocker arm and valve stem pin, the other end is rigidly fixed to 
ground. The spring is assumed to have a spring constant of 245 pounds 
per inch. The diameter is .07 inches; free length 46 inches. 
e. Valve 
This reciprocating valve is rigidly attached to the valve 
stem. The valve stem is assumed to be elastic. The upper end of 
the valve stem is connected to the common pin with the valve spring 
and right rocker arm. The diameter is .3 inches; length 3 inches. 
To complete the model two additional elements are added external to 
the system. Both provide for the Noxternal" clearances. They are: 

f. Cam 

The offset circular disk cam imparts sinosudial motion 
to the follower. 

g, Valve Seat 

This rigid seating surface for contact by the valve may 
be removed from the system when desired. 

The rocker arm, although modeled as two elements, is considered 
to be one continuous member. Thus only rotation about the fixed 
pivot point is permitted. 

2, Clearance Parameters 

The clearance or interfereance at the valve seat may be 
eliminated by the introduction of a predetermined parameter. Due tO 
the choice of the coordinate system, the parameter set less than the 


minimum value of the valve displacement has the effect of removing 
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the valve seat from the system. This provides for the Plescu bi laaby 
to conduct studies of the valve motion unimpeded by the valve seat. 

Additional flexibility of the model is provided by introduction 
of another clearance parameter. This parameter has the effect of 
eliminating the internal clearance petween the pushrod and left 
rocker arm. Thus these two members can be pinned together when 
desired. 

The circular disk cam rotates counterwise. This motion first 
lifts the follower by direct contact to its extreme position and 
next permits the follower to return to its initial position. The 
harmonic cycle is then repeated. For systems with cams of this type, 
the return motion of the follower must be accomplished by some force 
external to the cam. Therefore the valve spring is introduced to 
supply this force. For this work, simple harmonic motion is specified 
for the follower. The follower then oscillates about the midstroke 
position with an acceleration proportional to its distance from this 
position. The amplitude is half the stroke. The cam profile used 
to obtain the harmonic motion has been generated by the equation 
.5 sin (wt) inches, where » is the cam speed in radians per second. 
Similar studies of various other cam profiles are discussed in Refs. 


14-18. 


B, ELEMENT COORDINATES 

The model of the mechanism consists of two basic elements which 
account for all the relevant energy. The axial element, Figure Pape 
considers tension and compression forces. The bending element, Figure 


2.3, considers only pending energy. The coordinates of the five 
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Figure 2.2 AXIAL ELEMENT 
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Figure 2.3 BENDING ELEMENT 


elements of the system, i.e. pushrod, rocker arms, valve, and spring 


are chosen as indicated in Figure 2.4. 
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Figure 2.4 COMBINED ELEMENT COORDINATES 
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C. SYSTEM COORDINATES 

Assuming small deformation theory, the linear face! is composed 
of the above elements by superposition. The external members, the 
cam and valve seat, are also assigned coordinates for reference. 
As noted above, the pivot point and lower end of the spring are 
connected to ground and are therefore rigid. The systems coordinates 
fq} define motion. Since no motion is assumed for these two connections, 
no coordinates are necessary. Recalling that each member was assumed 
straight, the system coordinates will generally coincide with the 
element coordinates. The coordinate indices and their positive 


senses are shown in Figure 2.5. 
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Figure 2.5 FULL SYSTEM COORDINATES 


26 





D. DEFINITIONS OF CONFIGURATIONS 

The model shown in Figure 2.5 has two external clearances and 
one internal clearance. The term external refers to a displacement 
difference between a member of the system and an external member. 
Internal clearances are displacement differences between two internal 
members of the system. The clearances are: 

4. External To The System 

a. Clearance 1. The displacement difference between the 
cam, (coordinate nine), and the follower, (coordinate one). 

b. Clearance 3. The displacement difference between the 
valve, (coordinate seven), and the valve seat, (coordinate ten). 

2, Internal To The System 

a. Clearance 2. The displacement difference between the 
pushrod, (coordinate eight) and left end of the leit rocken sarm, 
(coordinate two). 

b. Clearance 4. The displacement difference between the 
valve stem, free end of the spring, and right end of the right rocker 
arm. The basic concept for the four clearance problem is presented 
in Appendix C. Since only the three clearance problem is considered 
for this work, clearance four is zero. Thus the members remain 
pinned and all move together at this connection. Only one coordinate 
is required to model the motion of the joint and is designated 
coordinate four. 

Thus the external clearances define positions where external 


forces must be considered, i.e. forces associated with coordinate 


ait 





one and coordinate seven, The internal clearances produce internal 
forces. These equal and opposite forces act upon the coordinate pair 
two and eight only when they are in contact. 

The three clearances, where positive displacement differences 
between the coordinate pairs may exist, will be referred to as noted 
above, i.e. clearances 1, 2, and 43. : Since clearance 4 will not be 
modeled, reference to this position will be referred to as coordinate 


four. Forces between clearances will be discussed in section 4.E. 


E. SYSTEM CONFIGURATION 
Depending upon the value of the force and displacement at these 
three clearances, the model can exist in one of eight possible 
configurations. If clearance 2 is considered pinned, the difference 
between coordinates two and eight is zero. The possible configurations 
is then reduced to four and is the study conducted by Anderson (Ref. 2). 
The eight configurations are diagramed in Table 2.1 and 2.2. 


They are: 





1. No contact at clearance 1, clearance 2, or clearance 3, 

2. Contact at clearance 33 no contact at clearance 1, or 
clearance 2. 

Be Contact at clearance 23 no contact at clearance 1, or 
clearance 3. 

4. Contact at clearance 2, and clearance 3; no contact at 
clearance 1. 

5. Contact at clearance 1; no contact at clearance 2, and 


clearance 4. 
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6. Contact at clearance 1, and clearance 3; no contact at 

elearance 2, | 

{. Contact at clearance 1, and clearance 23 no contact at 

clearance Ne 

8. Contact at clearance 1, clearance 2, and clearance 3, 

These eight configurations will be referred to as configurations I 
through VIII respectively. 

The physical constraints of the system require that at any given 
instant of time, one and only one configuration will be satisfied. 
The forces and displacement differences at clearances 1, 2, and 3 
are monitored to determine which configuration is valid for the 
calculated values. For example when the model is in a configuration 
that requires contact at a clearance and the calculated force at this 
clearance becomes zero (or negative), then contact is lost at that 
coordinate pair. Similarly, when the model is in a configuration that 
requires no contact at a clearance and the calculated displacement 
becomes zero (or negative), an interference exists and contact is 
regained at the coordinate. 

The configuration and the values monitored for the appropriate 
coordinates are listed in Tables 2.1 and 2.2. Since this model 
always has contact at coordinate four contact at clearance 2 will 
be ee to as contact at the "rocker", (between left rocker arm 


and pushrod). 
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F, SHIFTING OF CONFIGURATIONS 

As discussed above, the model shifts from one configuration to 
another, depending on first, the configuration the model is presently 
in; and second, the values of the monitored parameters. This shifting 
of configurations can be more easily understood by considering the 
flow diagram in Figure 2.6. The configuration schematic of the two 
clearance problem of Anderson (Ref. 2) appeared as a square, (face 3-4-8-7). 
As seen here the schematic of the three clearance problem becomes a cube. 
The solid lines and adjacent force parameter indicate a shift in the 
direction of the arrow caused by that particular force parameter. 
Similarly the dashed lines and adjacent coordinate parameters indicate 
a shift in the direction of the arrow caused by that particular 
negative displacement difference. Also indicated in the configuration 
flow diagram is an impact indicator. When a positive clearance becomes 
zero (or negative) a rapid change in the associated coordinate 
acceleration occurs which lends to high forces, i.e. impact. This 
is further discussed in section 3.F2. Impact is indicated by the 
rectangular boxes associated with the dashed shift lines. The figures 


within these boxes are the clearances at which the impact occurs, 
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Figure 2.6 FOUR CLEARANCE SYSTEM CONFIGURATION 
SHIFTING DIAGRAM 





IIl, MATHEMATICAL MODEL 


In this chapter the equations of motion which apply to the model 
are presented and analyzed. The development of the system mass, 
and stiffness matrices; the two eigenvalue problems; and the formulation 
of the system damping are discussed, An algorithm to shift the model 
to the configuration corresponding to the values of the monitored 
parameters is presented. Finally the numerical methods by which the 


solutions are obtained is explained. 


A. EQUATIONS OF MOTION 
The model is considered to be a linear viscous damping spring-mass 
system, excited by time varying forces. The differential equations of 


motion can be written in terms of the system coordinates as 


mi etae + 1c} {4} + [kK] fq} = {8 Gay 


where in the usual terminology, [M] represents the system mass matrix, 
[C] the system damping matrix, and [K] the system stiffness matrix. 
The formulation of these matrices has been developed using 
Lagrange's equations, with the matrix elements based on finite 
element techniques, Rubenstein and Hurty (Refs. 8 and 9). 
The elements of the model, as developed in section 2.B, have 
a discrete set of element coordinates (u, = 1,250.7 that represent 
the "n", rotational and translational element displacements. These 


element coordinates {u} are directly related to the system 


3h 





coordinates {q} by 


fu} =[B]{q} (3.2) 


The transformation matrix [8] is formed by applying the conditions 
of compatibility to the mechanism. For example, if element coordinate 
U4 of element three is rotated by one unit, recalling that linear 


theory is assumed, it is seen that the relation is 
abl = 
(uy)s = de (3.3) 


This scalar equation viewed as the matrix equation, (4.2), becomes 


(uy), =<000001000> | ?2 


Thus the fifth row of the transformation matrix [8] has been obtained. 
By similar procedures, the entire transformation matrix for the eight 
coordinate system is developed. Therefore equation (3.2) can now be 
written as equation (3.4). Where the column matrix of the element 
coordinates, {u}, are partitioned by their appropriate elements, 

and the 14 x 8 transformation matrix, [8] is partitioned into 
sub-matrices corresponding to the respective elements. A more 


concise form may be written for equation (3.4) as seen in equation 


.5). 
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(3.4) 
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B. MASS AND STIFFNESS MATRICES 


The element mass [M,] and stiffness [Ke] matrices are inherent 
characteristics of each element and can readily be formulated with 
finite element methods, (Refs. 1 and 8). The system mass [M) and 
stiffness [K] matrices can be formed with the aid of the transform 


matrix [8] as: 


D 
fh 
C00) -)) [8], (J, (89, (3.6) 
=1 


3, 
AL 
cK] -) > [8], (x,), (8), (3.7) 
i=] 


where fe), is the i*” element submatrix of [B); (KJ, is the i°® 


element stiffness matrix; and [Mo J; is the age element mass. The 


complete formulation is discussed in appendix A. 


C, DAMPING MATRIX 

The damping forces and hence system damping matrix is not as 
easily formulated since these values may depend upon the vibration 
of the system as well as elements exterior to it. To reduce the 
complexity of the solution of the differential equations of motion 
linear viscous damping is assumed. 

The system damping is formulated by Hurty and Rubinstein, 
(Ref. 9) by: (1) solving the eigenvalue problem associated with 


the equations of motion of the system in undamped free vibrations, 
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for the modal matrix; (2) using this modal matrix to form the 
diagonalized system mass and stiffness matrices; (3) from these 
matrices forming the diagonalized system damping matrix; and (4) using 
this diagonalized damping matrix to solve for the system damping 
matrix. 

The dynamic model consists of two eigenvalue problems. One 
eigenvalue problem is for the system of Figure 2.5. Here eight 
coordinates are required to solve the dynamic equations, since the 
pushrod and rocker arm are Eepeetny coordinate two and coordinate 
eight have independent motion. The other eigenvalue problem is 
associated with the model in contact at clearance 2 and takes the 
form derived by Anderson (2). When the motion at coordinates two 
and eight are equal, the eight system equations must be reduced to 
a lower order system by removal of the redundant equation to facilitate 
a solution. Thus the system damping matrix, as well as the system 
Mass and stiffness matrices, become order seven. This will be 
considered in section 3.H. 

1. Uncoupling the Equations of Motion 

When the differential equations of motion are derived from 
a set of convenient coordinates, the equations are normally coupled. 
If the system displacements {q} are expressed in terms of the normal 
modes, the governing equations will uncouple. Thus from the equations 


of free undamped vibration 


[M} {4} + [x] {a} = 0 (3.8) 
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the eigenvalue problem may be formulated as 


(twItu] - (k}) {a} = {0}. (3-9) 


The solution to this problem will yield the spectral matrix [2] 


and the modal matrix, [4]. 


il 
e 


cal (3219) 


e 


[3] 


£8hys (8}oe aoe (8), (3.11) 


The integer "n" is the number of degrees of freedom of the system 

i.e. n = 8 for the first eigenvalue problem discussed above and 

n= 7 for the second. Thel vj (i = 1, 2, ... n) are the natural 
frequencies (or eigenvalues) of the system; and the {3}. (4 =). $p2egeeee 
are the column vectors of mode shapes (or eigenvalues) corresponding 

to the natural frequencies. Now the governing coupled equation of 


motion for undamped forced vibration, 


(m] {4} + (K] {a} = {F3 (Geil?) 
may be decoupled by performing the coordinate transformation 


fa} = Le} in 


The dynamic equation can now be written in terms of the modal 


coordinates {Ni}, 


Raj ci} + Pel - 4 (3.13) 


SNS, 





the uncoupled, diagonal mass EN and stiffness Px] matrices, 


and the modal force vector {F *}, Hurty and Rubinstein, (Ref. 9). 


ee Relation of System Parameter To Modal Parameter 


The diagonal matrices, [mJ], aN and Pe] are each 


related to their corresponding system matrices by transformation 


Wath the modal matrix. 


Pa = C87" Or) £33 (3.14) 


Pe] = C83" (cl (81 (3.15) 
x J = [a7 (x) [81 (3.16) 


Additionally the modal acceleration {7} and displacement {1} are 
related to the system acceleration {q} and displacement {q} by the 


modal matrix. 


fq} = (3) (} (3.17) 


{qa} = [3] {7} (3.18) 


3- Calculation of the System Damping 


The system damping matrix may be obtained from equation 
(3.15) by matrix multiplication. Thus the diagonal elements of 
Pe J must be calculated before the desired result, [C], can be 
found. The governing equations are now uncoupled and the uncoupled 
mass Pad and stiffness Pk] matrices are known. Thus the solution 
formulated for the simgle degree of freedom system may be applied to 


the n degree of freedom system. 
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Equations (3.19) through (3.21), from Thomson (Ref. 3), may be 
utilized to obtain the desired diagonal elements of the uncoupled 
damping matrix Pe]. The natural frequency in radians per second 


of the undamped oscillation for a single degree freedom system is 


w= k/m (3.19) 


The damping factor (C), is defined as the ratio of the viscous 


damping (c) to the critical damping factor (c,). 
C= c/c, (3.20) 


The critical damping factor is also given by 


Cae 2d siti w = 2afm k (35.21) 


Upon substitution of equation (3.21) into equation (3.20), the 


viscous damping becomes 
c= 2CVknm (3.22) 
In terms of the uncoupled system, the value of the diagonal 


element of the natural frequencies, corresponding to equation (3.19), 


are 


m 
W, 11 
Ko9 
7s = wee (3-23) 
O ae g 
@ kon 
0 m 
Find 


where the k's and m's are the diagonal elements from lk] and Ee : 
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Similarly equation (3.22) can be written as 


C.. = 2c, k, 0, (1% 41, 2eeeer) (32225) 


es 


and thus the diagonal values of the uncoupled damping matrix Po] 
are obtained. 

An arbitrary choice of values for the damping factor vector {Cc} 
was considered on the basis of the natural frequencies. The damping 
factor, (S.), corresponding to the fundamental frequency is set 
equal to 4, The remaining elements of ic} — set to equal Wir 

Having established aN , the damping matrix [C] may be formed 
from equation (3.15) by matrix multiplication on the left by Ean 


and on the right by fay"! as 
eT IGT -) Get eh 5.25) 


This gives the desired result 


fo]- [eT Be Ey ae 


D. CONFIGURATION SHIFTING ALGORITHM 

For the two clearance problem four configurations were established. 
This work provided for three clearances having eight possible model 
configurations. For the addition of the clearance at coordinate 
four, a total of 4 clearances, there would exist sixteen possible 


system configurations. Thus the number of possible configurations 


he 





of an 'n" clearance system would be o,. SK generalized method of 
selecting the proper configuration on the basis of the monitored 
parameters is therefore necessary. 

The physical model will shift configurations only on the basis 
of two general parameters, displacement difference and force. Thus 
a fee, number system facilitates the problem of "computer book- 
keeping". This is accomplished by assigning the value one to a 
clearance that is in contact and a zero to a clearance not in contact. 
Por the unique choice of configurations established in chapter 2, 
Figures 2.4 and 2.5, the configuration matrix is developed as shown 


mefigure 3.1. 


CLEARANCES 


#1 #2 #3 





= |= |= |= OO 0 0 O 
—- =| © OF = =—& OC OO 
= O = 0 = O =| OO 


, 
2 
3 
4 
5 
6 
1 
8 


MOH PDA OMH?YAaAOOD 


Figure 3.1 Configuration Matrix 


The columns of this matrix are associated with clearances 
1, 2, and 3 as indicated. The rows are associated with the 
configuration number that the physical system has taken. Thus 
to retrieve the correct configuration the entire row is taken as a 


binary number and converted to a decimal number. Since it is not 
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convenient to have a configuration numbered zero, a one was 
arbitrarily added. 

For example, when no contact exists at all ees the 
binary number representing this configuration is 000. This binary 
number converts to decimal as 0. Upon adding the above mentioned 
value one, the decimal value becomes 1, which is the first row of 
the matrix and hence configuration I. 

The binary number 1 1 1 represents contact at all clearances 
and converts to decimal 7. Adding one yields the value 8, the 


th ; 
8°” row of the configuration matrix and hence configuration VIII. 





CONFIGURATION Vil CONFIGURATION V 
tep 1: ye 

Step 2: a” ja 

Step 3% 100 

Step 4: 100->V 


G) indicates binary addition 


Figure 3.2 Shifting Example No. 1 


The manner in which the program determines which configuration 
remmode? has shifted into is of primary concern. The developed 


algorthm is an extension of the above concepts, and best described 
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by way of examples. For the dynamic behavior pictured in Figure 3.2, 
the sequence of events is: 
(1) The model is in one of the eight given configurations, 
es. vol at time t. 
(2) The program calculates the required parameters for the 
next time step, ae Among these parameters are the monitored 
values. In this example the monitored value By has changed 
during the time step, 1.e. By has gone from Bo > 0 to F,< QO. 
(3) The negative force at coordinate two indicates the loss 


of contact at clearance 2. Thus the model now shifts to con- 


figuration V for time on (Figure 2.6). 


The four step shifting algorithm accomplishes this as outlined Tg 
Figure 3.2. The binary representation for configuration VII is 
retrieved from row seven in the configuration matrix, 110. Next 
the column of the binary number corresponding to the clearance that 
has been altered is removed. In this example, the element being 
the second digit of the binary number which correspond to the 
clearance 2 column. A one is then added, by binary addition, to 
this element. Carries are Fecegepsenaslel The new second digit O 

is placed into the binary number. Finally the new resulting binary 
number is converted to decimal and yields the new configuration 


number, V. 
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CONFIGURATION V CONFIGURATION VI 


Step 1: V-——>1 0 


SipeD 2: G 1 a 
Step 3: roy 
Step 4: 1.0 lee 


G indicates binary addition 


Figure 3.3 Shifting Example No. 2 


The second example for the shifting algorithm is shown in 
Pigure 3.3. The sequence of events for the dynamic behavior ‘is: 

(1) The model is in configuration V at time te 

(2) The program calculates the required parameters at the next 

time step ee In this example, the monitored displacement 

difference value a7 - do has changed from d7-F4 0 > 0 to 

17-440 <0; 

(3) The negative displacement difference at the coordinate 

pair 7 - 10 indicates that contact is established at clearance 

4. Thus the model now shifts to configuration VI for time 


ant 





The shifting algorithm accomplishes this as outlined in Keoure 5, 53 
The binary representation for configuration V is retrieved from row 
five in the configuration matrix, 100. Next the column of the 
binary number corresponding to the clearance that has been altered is 
removed. in this example, the element being the third digit of the 
binary number which corresponds to the clearance 3 column, The new 
digit, 1 is placed into the binary number. Finally the new resulting 
binary number is converted to decimal and yields the new configuration 


number, VI. 


BE, METHOD OF CALCULATIONS 

Newmark's B parameter method, (Ref. 4), forms the basis for 
the calculations. This one step numerical intergration procedure 
assumes the accelerations of the system coordinates are a linear 
function of time. Errors introduced by this assumption are given 
by wo A T/6, where ue is the highest natural frequency of interest. 
For small time steps this error is negligible. In this study the 


time step is taken as 
AT =i1x 107° seconds (527) 


It will be shown, by example that the solution to the equations of 
motion takes a convenient form. Precalculated matrix expressions 

May generally be used, resulting in large savings of computer 
execution time. A basic familiarity with the formulation of Newmark's 
method (Appendix B) will prove helpful and aid in the understanding 


of the following discussion. 
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The development of the problem calculations will assume the 
model has been operating in arSviciayspngonsttone contact at the cam 
and valve seat. The procedures are identical by ome to the 
remaining seven configurations. Thus the problem unknown quantities 


have been calculated for time t+ = t. The object is to solve the 


governing equations for the next time increment t = yntt 


1. Solution to the Governing Equations 
The full system matrices [M], [C], and [K] have been developed 
in sections 43.B and 43.C. The governing equations of motion are 
expressed as the matrix equation 4.1 in section 3.A. Written in 
concise form the governing equations for the eight system coordinates 


are 
[ms | (y,yor + fe; | ia + le, | fqsjtt! = (Ryo 
"(4 =4),2,..-83 j = 182,..58) (3.28) 


: and F a are unknown. The remaining 


ii 


= 2,3,..6,8) are equal to zero. Since 


; n+ 
Exterior forces, Fy 


applied forces ie ‘¢i 


the cam and follower are in contact, the known cam values at ref- 


n+1 . n+ n+1 


erence coordinate nine, dg 94g» and Io are equal to the 


oe N+ 
corresponding follower values at system coordinate one, q4 : 


—_ and ae Similiarly, the known values for the stationary 


. ee N+1 . ntl n+1 
valve seat at reference coordinate ten, dag 9 IO and dig are 


set equal to those corresponding values for the valve at system 


coordinate seven, e, ey and “me Displacements 
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mrt /. 
a. fi = 2, 3, -.-- 6, 8) are unknown. Expanding and partitioning 


equations 3.28 to take advantage of the known values yields 


ML 6 eee one 


BOG eee says 





n+1 n+1 


(3.29) 


©O 
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The next step creates two sets of equations. The first set, equations 
3.30 and 3.31 consists of individual equations in which the dis- 
placement, velocity ,and accelerations are known but the force is 

not. The second set, equation 3.32 consists of the matrix equations 

set equal to the zero force vector. In these equations all displacement, 
velocity ,and acceleration quantities are unknown. 


These equations of motion for the system are now written as? 
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The set of consistent system equations written in concise form are 


8 
F, = > (m, 4 ae + Cas al 1 kag ges (3555) 
it 
8 
F -)) (m5 cla + C75 ial " e, a) (3.36) 
i=1 
and 
2 
eo (N+ e nt n+1 nt) _ 
im J (aj P° + Cos.) CaP" + Cay] Cas} +) op oe = {0,3 (3-37) 
j=1 


Gis 2a. 6.0. Jee. bees O50) 
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The next step is the solution of the matrix equation 3.37 by 
Newmark's method, for the unknown accelerations. By integration the 


velocities and displacements are written in terms of the accelerations 


ca} = (a,Jo + Sao + Sa (3.38) 

and 
(a,}°*! = {a,}? + or (4,}" + * (a,}7 + aoe ie (3.39) 

(GIR 2) Dsus) 

Let 
fo}? = (4,9 + > (4) (3.40) 

and 
(8,}° = {a,}% + ar {a,}" + aCe (5.41) 


(gi = 2555 o- 6,8) 


Euostituting for { }” and {g}" equations 4.38 and 3.39 become: 


(a}t) = fo, J9 4 Saye (3.42) 
CA = (0, 0" t= sr be! su Gua) 


Gi - Qe oe O48) 


Next substituting for {4,]°t' ana {q,}"*1 into equations 3.37 and 


factoring yields: 


be 





2 : 
(Cm, J + to, = + ik, 4 == ) Al 


(3.44) 


2 
+ [c, 5] {o,}" 4+ [k, 5] {B,}° +) sean = {0,} 
Ise) 
(Giana 6,830 5 = 223, sa6,8) 
Let 
Ceerv, 1" = im, 5] - Lo, <4 a+ [kJ at (3.45) 


(lea) eo, e = 255, -- 6,8) 


Substituting for Caeav1~ in equation 4.44 and moving the knowns 


to the right side of the equation yields 


2 
Cmerv, 1 (4,}°*T = - Coz.) fa,}" - Ce, 5) 18,3" -) aa! | (5.46) 
k=1 


(i = 2,3,--6,8 3 j = Dan oolerey, 


Inverting the matrix Cemavy7! yields the expression for the unknown 


accelerations 
2 
(437 = Cao, We, oy)” - Camry, Ic, ,t8,)" - Ceav,,) () o,") 
k=1 
Ge= 2538666,983 5 =eeesee-6,8) (3.47) 


In equation 3.47 the non-time function matrix multiplications operations 


are precalculated as 


(pa, ,] = (ov, .] Ce, ,) (3.48) 
and 


[Dv 1 = [Eeiv, .] [k; 3] (3.49) 


(i = 2,3,..6,83 j = 2,3,-.6,8) 


Do 





The system accelerations, except the initial known values, are 


now calculated by 
2 : 
ee n+1 n+1 
ew 


(i = 2,3,226,8; j = ESOS), 


The system velocities and displacements are calculated from 


(a,}°*) = {o,} +S (ay (3.51) 
Z 
fa,}™*" = {8,} + (y,)") (3.52) 


(2 = 253556-6,8: j = 2,3,.6,8) 


The remaining system unknowns, Fy, and Fo may be directly 
memtcwlated from equations 3.30 or 3.35 and 3.31 or 3.36. This 
completes the calculations for the required solution to the problem 
for this time step. 

2. 'eeealculated Matrices 

Returning to equation 3.32 it is noted that the mass, damping, 
and stiffness matrices are unique to configuration VI. Similarly 
from equation 4.30 and 4.31 the mass, damping, and stiffness column 
vectors are also unique to configuration VI. Thus these matrices 
need only be calculated once if stored as layers to the full system 
Matrices. Figure 3.4 graphically indicates the matrix dimensions 
and elements for this and the remaining seven configurations. 
Additionally for the eight sets of precalculated matrices [M], [c], 
and [K] the matrix variables [BMIV], [DA], and [DV] given by equations 


5245, 3.48, and 4.49 may be precalculated for the standard AT and stored. 
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F, REFINEMENTS IN THE SOLUTION 
4. Incremental Time 
As seen in section 3.E, equation 43.27, the normal time step 
is established as 1 x 107° seconds. Generally the system configuration 
will change in a time step interval, rather than at either of the 
end points. When this occurs the present negative monitored value 


n+1 


calculated at ¢t t and the past stored monitored value that was 


calculated at t = t are used for a linear interpolation to obtain 
tne zero point, i.e. zero force or zero Tiguiaeemene difference. 
The AT, now designated as AT", corresponding to the zero point is 
used to calculate the system values based on the old configurations. 
If the system has more than one change in configurations for a given 
time step, the minimum AT' is used for the calculations, a ea. 
The newly calculated critical value is then monitored, and if still 
negative, the time is again reduced by a small amount until a positive 
value of the monitored parameter is obtained. Figure 3.5 summerizes 
the above procedures for a negative force parameter. 
2. Impact Forces 

Anderson, in Chapter IV (Ref. 2), puts forth a good account 
of the study of impact forces. Numerous references are cited. Thus 
the historical background of impact forces in relation to this study 
will not be discussed here. 

When the displacement difference values reduce to zero 
(or negative), contact at a clearance exists. This produces a 


step input to the velocity function and an impulse to the acceleration 
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function. Finite time steps are utilized to obtain solutions. 


Therefore one time step is sufficient to reach the step input value 


for the velocity curve, but two time intervals are required to 


account for the acceleration caused by the impulse, Figure 3.6. 
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Maus the acceleration at time t@*' is 


This accelation is readily calculated for impact at the cam or 
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IV. IBM-360 COMPUTER PROGRAM 


This computer program is coded in FORTRAN IV, (Refs. 10-13), 
to run on an IBM-360/67, with 0S/360 release 18. The simulation 
is modeled after basic concepts put forth by Anderson, (Ref. 2). 
The original model program has been completely restructured anc 
formulated to yield general procedures for the analysis of the | 


dynamics of an elastic link mechanism with an internal clearance. 


A. GENERAL DESCRIPTION 

The program requires 234 K bytes of core storage for 850 degrees 
of cam rotation, including the CALCOMP graph buffer. The offline 
printer array (GRID) is not included in this core storage figure. 
Execution time is a function of the speed of rotation of the cam. 
For 11,000 RPM with 850 degrees, execution time is about 4.43 minutes. 
The primary output parameters fill arrays for SYSOUT = B and 
sYsOULT = C. The B queue is for transfering data, via punched cards, 
for use in the SDS-9300 computer. The C queue is for plotting 


graphs via the CALCOMP. 


B, GENERALIZATION MATRICES 

As noted in section 3.D, the number of system configurations is 
, where "n" is the number of clearances. Therefore, the most 
general computer programming was utilized, where possible. Several 
arrays of prestored computer "instructions" are either poner 


generated or formed by the user. These arrays set looping variables, 
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dynamically dimension work arrays, and establish program calculation 
flow indicators. 
{. Matrix CONFIG 
The 8 by 3 matrix CONFIG, Figure 4.1 was described in section 

4.D. This integer array indicates the physical configurations that 
the model has oe based on the value of the monitored parameters. 
The matrix is formed by the program. The user must supply the value 
for the maximum number of configurations the system is permitted to 


take. 





Figure 4.1 Matrix CONFIG 
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2. Matrix IDM 

The 8 by 1 column vector IDM indicates the dimension of the 
Matrices that are used for calculations. This array is naGlevecieerl 
in Figure 4.2 for the three clearance system. The program parameter 
ISYS is the matrix dimension of the equivalent system mass, stiffness, 
and damping matrices and for the calculated program matrix variables 
BMIV, DA, and DV. The dimensions of the various layers correspond to 
their configuration number and to the rows of IDIM, i.e. when the 
system is in configuration I, layer one matrices meapa dae dimension 
indicated in row one of JDIM. This integer array is calculated by 


Ld 


the program. 


ee) 


NON ns ~~ 


ON 


Figure 4.2 Matrix IDIM 


Db. MMetrix WAYER 
The 8 by 1 array LAYER, Figure 4.3, is an integer column vector 
supplied by the user. The elements of the vector consist of ones 
or threes. As the name implies, the layer of the mass, stiffness, 
damping, BMIV, DA, and DV matrices used in the calculations correspond 


to the configuration number the system has taken. The eigenvalue 
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problem for contact at the internal clearance (layer 1) or no 

contact (layer 3), determines the elements of LAYER, Configurations 
TiI, IV, VII , or VIII for contact at the rocker—pushrod, use the 
same eigenvalue problem for the reduced system, i.e. the damping 
matrix is calculated with layer 4 matrices. Configurations I, II, V, 
or VI for no contact at the rocker-—pushrod, use the full system 
eigenvalue problem. The damping matrix is calculated using layer 1 
matrices. Figure 4.3 shows these values of LAYER for the three 


clearance model with two eigenvalue problems, 


NNN eS ON ON 


Figure 4.4 Matrix LAYER 


4. Matrix IQ 
The 8 by 8 integer array IQ, Figure 4.4 indicates the system 
coordinates for which the Newmark routine will calculate the unknown 
quantities i.e. displacements, velocities, accelerations, and forces 
for a given time step. The columns of the array correspond to the 
configuration number except column 1 and 4. These are additional 
columns used with configurations VI and VIII respectively. The 


array 1S user supplied. 
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foe Matrix LDP 
ims Dy 5 integer array, Figure 4.5, consists of zero 
elements except where impact may be expected to occur at the externas 
clearances. Here again the rows correspond to the configuration 
number and the columns correspond to the clearance number. The 
Matrix IDP is used in the calculation of AP as a portion of Newmark!s 


method. The array is calculated by the program. 
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Figure 4.5 Matrix IDP 
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bee Watrix ICORD 

The 3 by 2 integer array ICORD, Figure 4.6 indicates the 
coordinate pairs of the internal and external clearances. The 
elements are stored row wise by pairs. The rows correspond to the 
clearance numbers 1, 2, and 43 respectively. The first colum 
corresponds to the index of the system coordinate that has the greater 
positive displacement when no contact exists at the clearance in 
question. The second column is assigned to the remaining index of 
the clearance pair. The array ICORD is user generated, Figure 4.6 


shows the matrix for the three clearance problem defined in chapter 2. 


1 2 
2 8 
T 10 


Figure 4.6 Matrix ICORD 


7, Matrix IMPACT 
The 3 by 1 column vector IMPACT is program generated. ‘The 
integer elements are continuously modified during program execution. 
The rows correspond to the clearance number pairs. Values 1, 2, 
or 7, of the system coordinates are stored when impact occurs for 


a given time step at a corresponding clearance. 


C. DYNAMIC DIMENSIONING 

The principle of dynamic dimensioning is widely used in advanced 
programming and often misunderstood. Thus a short example of its 
concepts as applied to this work is presented here. Reference to Figures 


4.7 and 4.8 will aid in the understanding of the concepts. 
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Figure 4.7 Dynamic Dimensioning of Array SYSM 
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ARRAY BMIV 


subroutine: THES IS BKS TEP DEL 


DIMENSION: (8,8,9) (8,8) (8,1) 


CALI PRE TEP (BMEV (edi 9 ) 7 


CALL DEL [BMIV,...] 





Figure 4.8 Dynamic Dimensioning of Array BMIV 
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In subroutine THESIS the mass array SYSM and the calculation 
array BMIV are exactly dimensioned indicating the total number of 
elements which are stored by columns in the computer. When eet ine 
BKSTEP is called, the entire array SYSM is passed. However the 
array BMIV is initialized in the call as the first value of the 


last layer, in this example. This is shown as: 


Sueroucine THESIS 
DIMENSION SYSM (8,8,8), BMIV(8,8,9) 


CALL BKSTEP (SYSM, BMIV(1,1,9), -.-) 


Thus in subroutine BKSTEP, the array SYSM must be dimensioned with 
the exact number of rows and columns. Since the computer has 
already provided for the total space of the array, the third 
dimension is dynamic and if entered may be any positive number equal 
to or less than the originally assigned dimension. Since the entire 
array 1s utilized the third dimension is retained. For the array 
BMIV only the last layer is used, thus the third dimension is 
dropped. 

Subroutine BKSTEP 

DIMENSION SYSM(8,8,1), BMIV(8,8) 


CALL DEL (SYSM(1,1,C), BMIV, ...) 


When subroutine DEL is called from BKSTEP the array, SYSM is 
now initialized with the variable C, (C = 1,2,...8). BMIV, now only 
considered an 8 by 8 array, is not initialized here, The dimensioning 


in subroutine DEL is accomplished similar to that done in BKSTEP. 
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The array SYSM must again be dimensioned as above with the third 
dimensioning parameter being dynamic. The array BMIV is used as a 
double dimensioned array in subroutine DEL. By arguments similar 

to those above, the columns of the array may now also be dynamically 


dimensioned as: 


subroutine DEL 


DIMENSION SYSM(8,8), BMIV(8,1) 


Figure 4.7 for array SYSM and Figure 4.8 for array BMIV graphically 
indicate the manner in which these arrays are passed from the 


subroutines, THESIS—BKSTEP—DEL, 


D. PROGRAM STRUCTURE 

The main program only sets up the plotting arrays and then calls 
subroutine THESIS, the main routine that controls the program execution. 
The calling order is presented in Figure 4.9. The CALCOMP graphing 
subroutines indicated by an asterisk are supplied by the computer 


center. 


BK. SUBROUTINES 

The function of each of the subroutines and a description of 
their execution are included in this section. The concept of 
execution of object modules stored on disk was used to a large extent. 
This accounts for the large number of subroitines and their limited 
content. Flow graphs for the majority of the routines are collected 
and presented at the end of this section. Hach Pateatine is listed 


in its entirety in the IBM-360 Computer Program Listing chapter. 
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1. Subroutine THESIS 

This subroutine controls program execution, establishes 
calculation arrays (MCK, CONMAT, MATRIX, DEL), initialing parameters 
(INPUT, INITAL) and contains the major time step loop, AT = 107° 
seconds. As seen in Figure Leo}: most subroutines are called from 
THESIS. Figure 4.10 indicates the flow diagram. 

The time parameter wt advances by an amount wAT. The previous 
values of displacement, velocity, acceleration, and force are stored 
(ADJUST). The new cam values are calculated (CAM) and the know 
system parameter are initialed according to the present configuration 
(BRANCH). Newmark's method (NEWMARK) is then employed to calculate 
the remaining unknown system quantities. The system configuration 
value is next stored prior to checking the newly calculated quantities 
(MONTOR) for any change in the configuration number. 

The system is now checked to see whether an impact has occurred 
and at what coordinate index. If the configuration has changed, the 
quantities are recalculated (BKSTEP) for the time that the configuration 
changed. Then force values are reset (FIXFOR) to coorespond to the 
system configuration. 

The amount of degrees of cam rotation is calculated. If the 
next whole degree has been reached the system values are output for 
use as the graphics (AGT) and draw (ARRAY) data. Next if the required 
number of degrees have been evaluated the desired graphs are plotted 
(GRAPH). Otherwise, the program branches to the start of the time 


step loop (statement number 360) and the cycle is repeated. 
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2. Subroutine BKSTEP 
This subroutine is similar to subroutine THESIS in that all 
program variables are recalculated to correspond to the time at which 
the system configuration changed. Once this routine is entered only 
the parameter corresponding to the configuration shift is monitored 
until a positive value is calculated. Figure 4.11 indicates the flow 
diagram. 
3. Subroutine MCK, DVVS, and DINVS 
This subroutine originally coded by Anderson and Winfrey, 
was modified to account for the additional eigenvalue problem. The 
sizes and material properties of the system members are assigned. 
The element of the full and reduced system mass and stiffness are 
assigned to their respective arrays, SYSM and STF. The damping 
ratio, (ZTA) is set on the basis of the fundamental frequency of 
the system being considered. The eigenvalue problems are solved by 
subroutine DVVS. Subroutine DINVS is called to obtain matrix inverses.. 
4. Newmark's Method Subroutines 
Subroutines NEWMRK and DEL are the primary program calculation 
routines, and considered together comprise Newmark's method. The flow 
diagram for subroutine DEL is seen in Figure 4.12 and for subroutine 
NEWMRK, Figure 4.13. 
(a) Subroutine DEL 
Subroutine DEL calculates the eight layers of the matrix 
arrays BMIV, DA, and DV, based on the full time step AT. These 
matrices are precalculated and stored. When DEL is called from 
subroutine BKSTEP, an additional layer of the above matrices is calculated 


as required, with a partial time step, AT sn Palme seconds. 
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(b) Subroutine NEWMRK 
Subroutine NEWMRK calculates the coordinate accelerations 
by Newmark's method. The theory of this method is presented are 
appendix B. Its relation to this study is presented in section 3.E. 
After the system accelerations are known, the displacements and 
velocities can be calculated. Finally the forces are calculated at 
clearance coordinates 1, 2, and 7. 
5. Subroutine for Shifting Configurations 
The two subroutines that provide the program the ability to 
shift from one configuration to another, based on the values of the 
monitored parameters, are MONTOR and SHIFT. Figure 4.14 indicates the 
flow diagram for subroutine MONTOR. The subroutine that determines the 
amount by which the time step should be reduced, after a change in the 
configuration has been indicated, is DIPRIM. The flow diagram for this 
subroutine is seen in Figure 4.16. 
(a) Subroutine MONTOR 
Subroutine MONTOR checks the forces and/or the displacement 
differences at the three clearances to determine if any has become 
negative. As indicated by the shifting flow diagram, Figure 2.6, only 
three quantities must be monitored at each clearance, either forces 
and/or displacements differences. When at least one of the values 
becomes negative, subroutine DTPRIM is called. 
(b) Subroutine DTPRIM 
subroutine DITPRIM linearly interpolates for a best 
estimate of the time the monitored value became negative, i.e. when 


the system changed to a different configuration. The calculated time 
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is stored in array DTP. If more than one change to the system 
configuration occurs in a single time step, the first change is 
chosen, i.e. the minimum value among the elements of the array DTP. 
The clearance index number that causes the system to change is then 
sent to subroutine SHIFT. 

(c) Subroutine SHIFT 

Subroutine SHIFT, described in section 3.D returns the 
new configuration based on the parameters sent by MONTOR. 
6. Subroutine INPUT 

This subroutine reads as data the formal parameters that are 
frequently varied. Figure 4.17 presents the flow diagram for this 
subroutine. These parameters include revolutions per minute (REM) 
of the cam, clearance indicators (CLEAR, PIN 2&8) printing indicators 
(MPRINT, WIDEG, KPRINT, NCOUNT, IPRINT), output data indicators 
(IDRAW, AGTOUT, CARDS, TAPE IPLOT) and the tipe step (DELT). Graph 
titles are read into the arrays XNAME, YNAME, LABEL, and TITLE. 
The array TITLE is modified to indicate programs without a valve 
seat. The array LABEL is modified to indicate a two clearance problem. 
The lengths of the member elements are assigned (REALTH, AGTLTH) and 
the lever arm of the rocker (ROCKER) is calculated. 

7. Subroutine CONMAT 

This subroutine establishes two of the program generalization 
arrays, CONFIG and IDP. Figure 4.18 shows the flow diagram. The 
routine is called only once during the program execution. The coding 
for both arrays may be easily expanded to provide for problems with 


more than three clearances. 
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8. Subroutine INITAL 

This subroutine initiates the starting values of displacement, 
velocity, acceleration, and force for the system coordinates. The 
starting configuration VII, requires that layer seven matrices be used 
for all calculations. Subroutine DINVS is called to establish the 
inverse stiffness matrix necessary for the displacement calculations. 
The program variables for the follower are set equal to the cam values 
and those for the left rocker arm are set equal to the pushrod values 
by calling subroutine ADJUST. 

9. Subroutine CAM 

This short subroutine calculates the displacement, velocity, 
and acceleration at the cam (coordinate 9) for the given wt. The flow 
diagram is presented as Figure 4.15. The displacement is given by 
the equation: 

Io = .5 sin wt + XZERO 
This routine is called at least once, every time increment from 
subroutine THESIS. Additionally cam is called by subroutine BKSTEP 
for intermediate values of wt. 
10. Subroutine BRANCH 

This subroutine is called prior to making calculations for 
a new time step. Figure 4.19 indicates the flow diagram. The eight 
sections of the routine initialize the system parameters corresponding 
to its present configuration. Additionally the dimension parameter 
ISYS is set from the array IDIM. A correction ot the above established 


acceleration is accomplished when the array IMPACT indicates that a 
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clearance displacement difference is zero. The adjustment of the 
acceleration produces a corresponding increase in the later calculated 
merce. 
11. Subroutine ADJUST 
This short subroutine equivalences the displacement, velocity, 
and acceleration quantities at a system coordinate to those at another 
system coordinate. For example, when contact exists between cam and 
follower the above named quantities for coordinate one are set equal 
to those for coordinate nine. Figure 4.20 presents the flow diagram. 
12. Subroutine FIXFOR 
This subroutine is primarily used to zero the calculated 
forces on the basis of the system configuration. This is necessary 
to eliminate small negative or small positive forces from being 
plotted. This routine has no effect on the program calculations. 
Its flow diagram is seen as Figure 4.21. 
13. Output Subroutines 
(a) Subroutine AGT 
The subroutine AGT outputs data on punched cards or 
magnetic tape for use in the SDS-9300 program calculations of the 
graphic display. The flow diagram is indicated in Figure 4.22. This 
Subroutine first adjusts the system displacement coordinates by an 
amount XZERO so as to present the initial position of the graphics 
display as horizontal. In addition to the displacement coordinates, 
reference values are sent for the cam rotation and the valve seat 


position. 
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(b) Subroutines ARRAY, GRAPH, GPRINT, OLPLOT 
The subroutine ARRAY is called for the first value of 
each whole degree and establishes the array FPLOT, which is used LOT 
the plotting of graphs on the CALCOMP or offline printer. Figure 4.23 
presents the flow diagram. Subroutine GRAPH is the last subroutine 
called by THESIS after the desired number of degrees has been evaluated. 
This subroutine, as seen in the flow diagram, Figure 4.24, organizes 
the calls to the offline printer subroutine OLPLOT for IPLOT equal to 
YES and/or the CALCOMP subroutine GPRINT for IDRAW equal to YES. The 
flow diagram for subroutine GPRINT is indicated in Figure 4.25. 
Subroutine OLPLOT is a modified version of the Computer Center crovided 
routine UTPLOT. One of the modifications provides for an expansion 
of the graph output to a maximum of six pages. The number of rezes 
desired is assigned to the program parameter, IP, by the user. 
Additionally the users must adjust the dimensions of the main crogram 
array, GRID to conform with the desired pages. GRID must be dizensioned 
as LINES by 81 where the dimension parameter LINES is: 
LINES = (80)(IP) +1 
(c) Subroutine DRU-- 
The many subroutines beginning with DRU-- were written 
for plotting graphs on the CAICOMP plotter by Winfrey. This package 
also contains computer center provided routines, indicated by en 
asterick in Figure 4.9. 
(d) Subroutine LIST (DEBUG) 
Subroutine LIST outputs the numerical formal parareters 
on the offline printer for either "debugging" purposes or as reference 


values when viewing the graphs/graphics. 
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(e) Subroutine PRNTMX 
Subroutine PRNIMX outputs all the system matrices, 


when the program parameter MPRINT is equal to YES. 
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Figure 4.10 Subroutine THESIS Flow Diagram 
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Figure 4.11 Subroutine BKSTEP Flow Diagram 
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Figure 4.12 Subroutine DEL Flow Diagram 
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Figure 4,13-1 Subroutine NEWMRK Flow Diagram 
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Figure 4.14 Subroutine MONTOR Flow Diagram 
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Figure 4.15 Subroutine CAM Flow Diagram 
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Figure 4.16 Subroutine DIPRIM Flow Diagram 





Subroutine 
INPUT 


LABEL 
TITLE 

REM 

CLEAR, PIN2$8 
MPRINT 


YNAME. Gis" eee. 31 OD 
XNAME 





WRITE 
INPUT VALUES 
IS WRITE 
OUTPUT ON . BOF 
TAPE REWIND 


Figure 4.17 Subroutine INPUT Flow Diagram 
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Figure 4.18 Subroutine CONMAT Flow Diagram 
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Figure 4.19 Subroutine BRANCH Flow Diagram 
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Figure 4.20 Subroutine ADJUST Flow Diagram 
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Figure 4.21 Subroutine FIXFOR Flow Diagram 
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Calculate array QR4 values 


for graphic display 








Funch card Write tape 


output = | output 


Figure 4,22 Subroutine AGT Flow Diagram 
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for graph output 





Figure 4,24 Subroutine ARRAY Flow Diagram 
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Figure 4.24 Subroutine GRAPH Flow Diagram 
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Figure 4.25 Subroutine GPRINT Flow Diagram 
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V. SDS-9300 COMPUTER PROGRAM 


This computer program is coded in SDS FORTRAN IV, (Ref. 5-1), 
to run on a Scientific Data System Model 9400 computer interfaced 
with an ADAGE Graphics Display Terminal, Model 10 (AGT-10). This 
program provides the user the capability to observe the dynamic 
response of 360 degrees of cam rotation of the awe train assembly. 
The displacement data generated by the previously described program 
is the major input. Finite element techniques are used to reformulate 
the input displacement data to provide for the pictoral responses 
of bending and axial deformations. Due to the limited core storage 
of this computer, overlaying techniques have been employed. | 
A, METHOD OF ELASTICITY CALCULATIONS 

As shown in Figure 2.4,the system consists of five elements. 
Three of these elements, (1,4, and 5) are axial members. The 
remaining two, (2 and 3) are bending members. Thus the elastic 
mode shapes for each of the elements as a function of the distance 
along the member will yield the required deformations. R, C. Winfrey, 
reference (1) Chapter V, has done considerable work in the formulation 
of the required mode shapes by finite element techniques. The mode 


shapes for the axial elements as a function of x are: 


8, (x) (5.1) 


2 &) (5.2) 


where L is the longitudinal dimension of the bar. 


o> 





The mode shapes for the bending elements as a function of x are: 


#, @) = 5 @? - beh + 2°) Ce) 
tp) = GAL - 20°) (5.4) 
#3 (x) = GaP - 2x8 n+ x) (5-4) 
#4 (x) = 5 (EL #) (5.6) 


where L is the longitudinal dimension of the bar. 

The system displacements are calculated by the IBM-460 program 
and sent as data to the SDS-9300. These displacements and the above 
mode shapes are superimposed to yield the required rigid body motion 


and elastic deformations. 


B, PROGRAM STRUCTURE 


1. General Description 





The main program controls the program execution. The calling 
order is presented in Pigure 5.1. The subroutines indicated by 
an asterick are supplied by the computer laboratory, and are used to 
Manipulate and initialize data sent to the ADAGE Graphics Terminal. 
The SDS computer total core available for program execution, 
exclusive of the control and real time monitoring functions, is 
approximately 18 K bytes. Since the total program length greatly 
exceeds the core available, a program overlay technique consisting of 
storing only active subroutines during execution is employed. This 
achieves maximum use of the main core storage and allows the program to 


execute. 
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Figure 5.1 
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The overlay tree structure is depicted in Figure 5.2. Although 


necessary for this program, overlaying significately increases the 


program execution time. 


MAIN READ BARL BAR2 BAR5 BAR4 


ZZ 
INITAL SPROTR MIDBAR DiSBEY JOINT 


1. Remains in core 
2. Overlay structure 


Figure 5.2 Overlay Structure 


The bowing and thinning of the axial elements, as seen on the 

_ graphics display photographs, have been added as an aid in observing 
the deformations of these members. The mode shapes given by equations 
5.1 and 5.2 are only valid for the elastic line of a one dimensional 
bar. However the small uniform extension and compression of the 
element were impossible to satisfactorily reproduce on the display. 
Thus the bowing and thinning were added as a function of the axial 
deformation. The thinrer the bar, the more axial compression it 

has undergone. Conversely, the thicker the bar the more axial 


extension it has undergone. 
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2. Sense Lights 


The SDS computer being a hands-on system, has several physical 
sense lights. When lit, logic within the program may alter Bete 
instructions as desired. 

a sense Light 1. 

Activation of sense light one branches the execution to 
the head of the program. Hxecution then continues and variables are 
reinitialized as coded in the program. New data points may also be 
read in, displaying a different set of displacement data. 

b. Sense Light 2. 

During program execution this sense light indicates that 
the display will utilize ADAGE Graphics Terminal, number two. When 
not lit, the SDS computer is interfaced with the AGT unit number one. 

c. sense Light 4. 

This sense light branches the program to a convenient 
point to reinitialize the originally coded program parameters. This 
is accomplished by activation of the SDS control console teletype 
from which the user may input any program variables that appear in 
the namelist statement. 

d. Sense Light 6. 

Activation of sense light sixcalls the subroutine TEXTR 
for each cam rotation increment. This enables the motion of the 
display to be halted for further study at any point of interest in 
cam rotation. Depression of the return key on the AGT terminal will 


advance the display to the next display of cam rotation, 


Ne, 





C. MAIN PROGRAM 

The main program initializes parameters, establishes the data 
arrays and contains the calculation loop. Prior to the execution 
of the program the desired sense lights should be activated. The 
program is halted almost immediately aiter execution by an INPUT 
statement. This statement allows the user to have positive control 
over assignment of the AGT unit. The designation of the wrong 
unit at this point will case a core image failure. The program 
then must be reloaded, Therefore at this point the 9300 teletype 
Control console should be checked to insure the correct unit is 
listed in the message to the user. To proceed with the program 
execution a "*" must be typed on the control console followed by 
the carriage return key. 

Subroutines are called to display the fixed elements of the 
model. The two arrays that reproduce the circular cam are pre- 
calculated. Next the frequency ratio is calculated and displayed. 
The cycle loop is then entered where subroutines are called to 
display the movable elements. This loop is incremented by a 
specified amount of cam rotation until 360 degrees of rotation 
has been reached. At this point the cycle is repeated. For each 


cam increment, the sense lights are tested. 


D. SUBROUTINES 
The function of each ot the subroutines and a description of 
their execution are included. Hach subroutine is listed in its 


entirety in the SDS-94300 Computer Program Listing. 
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1. subroutine INITAL 
This subroutine initializes the text and graphics directory 
arrays, DIINIT and DGINIT, respectively. The text response, TEXTR 
is initialized to block number one at display location 39, 46. The 
title display is encoded and presented on the AGT. The graphic data 
blocks are initialized to zero, assigned the appropriate tirst word 
header (IHEAD), and then assigned sequence numbers in DGINIT for 
later use. 
2. subroutine READ 
This subroutine reads the title and label for the program 
and prints them on the ofttline printer. The program parameters 
are read, (XZERO, AMP, RPM, WN, ISEAT, and PIN28) and also printed 
for reterence. The user is then questioned trom the AGT display, 


Mmoto 5.1. 






LISTING OF COORDINATES DESIRED 
COMMAND: YES _ 


Photo 5.1 AGT QUESTION 


An atfirmative response will list the IBM-360 displacement data on 
the otfline printer. The response is entered via the graphics 
terminal. 

The displacement data is next read from magnetic tape or 
punched cards, depending on the control specification at the time 
the program was loaded. The array parameters read are, Q1i through 


Q9, OMEGA, and CONFIG. The Q's are as defined in the IBM program but 


nO 





modified by the amount XZERO, OMEGA is the amount of revolution of 
the cam in radians, CONFIG is the corresponding configuration 
number. Next the displacements are printed, if desired. Finally 
the title screen is nulled and control is returned to the main 
program. 
Se ocubroutine SPRCTR 
This subroutine is only called once for a given problem. 
The display for the rigid portion of the spring is calculated and 
displayed. The valve — position and geometry is calculated. 
The parameter ISHAT is then checked to determine whether the problem 
considers contact at the valve seat. For those problems with no 
contact, the valve seat location is adjusted to better picture the 
position. Finally the valve seat is displayed and control returns 
to the main program. 
4. Subroutine MIDBAR 
This subroutine is called once. Subroutine MIDBAR calculates 
and displays the fixed pivot point. 
5. Subroutine DISPLY 
This subroutine is called once during each cam increment. 
The function of this routine is to display the text for the system 
Configurations. 
6. Subroutine CAM 
subroutine CAM is called once during each cam increment. 
This routine calculates the cam position, given the cam rotation 


(OMEGA ) and displacement of the top of the cam (q9)- 


102 





7. Subroutine BAR 
Mis Se ae is called once during each cam increment. 
Seerdinates q., and dg are used to calculate the position of the 
pushrod. Equations 5.1 and 5.2 calculate the deformations. The 
pushrod is then displayed on the AGT. 
8. Subroutine BAR2 
This subroutine is called once during each ‘cam increment. 
Coordinates dos Iz and des and equations 5.3, 5.4, and 5.5 are 
used in the calculations of position, oo wer and deformations 
of the lett rocker arm. The parameter PIN28 is checked to determine 
whether the pin should be displayed at clearance two. The left 
rocker arm and pin, if required, are displayed on the AGT. 
9. Subroutine BARZ 
This subroutine is called once during each cam increment. 
Coordinates Vy Qs» and des and equations 5.3, 5.4, and 5.5 are used 
in the calculation ot position, geometry, and deformations of the 
right rocker arm. Calculation of the position and image of the 
rocker arm- valve stem- valve spring joint pin is accomplished by 
subroutine JOINT. The right rocker arm and pin are then displayed 
on the AGT. 
10. Subroutine BAR4 
This subroutine is called once during each cam increment. 
Coordinates q) and Qo are used to calculate the position ot the valve 
stem. Equations 5.1 and 5.2 calculate the deformations. The valve- 


valve stem is then displayed on the AGT. 
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11. Subroutine JOINT 
This subroutine is called once for the right rocker arm pin 
ema, if required, once for the left rocker arm pin, during each cam 
increment, The pins! position and image are calculated and displayed 


on the AGT, 
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A. MODE SHAPES AND NATURAL FREQUENCIES 

The mode shapes and natural frequencies presented in Figure 6.1 
result from the eigenvalue problem for the full eight coordinate 
system, Figure 6.2 shows the values resulting from the seven 
coordinate reduced system, for contact at the rocker. 


NATURAL FREQUENCIES w, (EIGENVALUES) x 10° 


Figure 6.2 REDUCED SYSTEM MODE SHAPES AND NATURAL FREQ. 
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B. PRE-CALCULATED SYSTEM MATRICES 

The remaining system matrices are presented and grouped according 
to their layer (or configuration) number. The mass, stiffness, and 
damping matrices for layer one are unaltered and directly correspond 
to the full eight coordinate system eigenvalue problem. The mass, 
stiffness, and damping matrices in layer three are unaltered and 
directly correspond to the reduced seven coordinate system eigenvalue 
problen. 

The pre-calculated matrices COLM, COLC, and COLK located in layers 
one and three are associated with the calculations for configurations 
six and eight respectively. Thus, these matrices are not listed in 
their correct program layer, but with the correct configuration with 
which they are used. 

The order of presentation within each layer is: (1) Mass, 

(2) Damping, (3) Stiffness, (4) COLM, (5) COLC, (6) COLK, (7) BMIV, 


(8) DA, and (9) DV. 
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Figure 6.3 CONFIGURATION I - Mass, Damping, and Stitiness 
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Figure 6.4 CONFIGURATION I - Precalculated Matrices; BMIV, DA, and DV 
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Figure 6.5 CONFIGURATION II - Mass, Damping, and Stiffness 
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Figure 6.7 CONFIGURATION II - Precalculated Matrices; BMIV, DA, and DV 
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Figure 6.8 CONFIGURATION III - Mass, Damping, and Stiffness 
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Figure 6.9 CONFIGURATION III -— Precalculated Matrices; BMIV, DA, and DV 
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Figure 6.10 CONFIGURATION IV - Mass, Damping, and Stiffness 
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Figure 6.12 CONFIGURATION IV -— Precalculated Matrices; BMIV, DA, and DV 
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Figure 6.13 CONFIGURATION V - Mass, Damping, and Stiffness 
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Figure 6.15 CONFIGURATION V — Precalculated Matrics; BMIV, DA, and DV 
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Figure 6.16 CONFIGURATION VI - Mass, Damping, and Stiffness 
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Figure 6.19 CONFIGURATION VI - Precalculated Matrices; BMIV, DA, and DV 
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Figure 6.20 CONFIGURATION VII - Mass, Damping, and Stiffness 
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Figure 6.22 CONFIGURATION VII - Precalculated Matrices; BMIV, DA, and DV 
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Figure 6.23 CONFIGURATION VIII ~ Mass, Damping, and Stiffness 
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Figure 6.25 CONFIGURATION VIII - Precalculated Matrices; BMIV, DA, and DV 





C. SAMPLE PROBLEMS 

One set of results for 11,000 RPM was obtained in CALCOMP graph 
form and as AGT graphic display photographs. At this high speed 
of cam rotation, the pushrod is frequently thrown off the cam. “he 
results are completely analyized with numerious reference points 
plotted on the graphs. 

The second example problem was chosen at a slower speed, 9,000 
RPM. At this speed the pushrod remains in contact with the cam. 
A third intermediate example, 10,000 RPM was chosen for a comparison 
between the two above examples. Figure 6.26 shows the relationship 
between the sample problems. 

Photographs from the AGT display are presented after the CALCOMP 
graphs for each case of example one. It is noted that large scale 
factors are introduced for both bending and axial deformation as an 


aid to the user. 


THREE TWO 
CLEARANCES 





Figure 6.26 SAMPLE NUMERICAL RESULTS 





1. Example 1A, Three Clearance, 11,000 RPM, Valve Seat 


This analysis will commence at 270 degrees of cam rotation. 
One cycle will be discussed, from point (A) to point (AA). The 
following figures are utilized for the analysis: 

a. Figure 6.27 


(1) Displacement difference between cam and pushrod 


=? 
(q, = Ig) x 10 


(2) Displacement of lower end of pushrod i.e. follower, 


— 
q, x 10 


(3) Cam Profile Gg x 107" 
be emuemie 6.28 


(1) Displacement difference between pushrod and rocker 
-2 

(a, a Qe) ea 
(2) Displacement difference between cam and pushrod 


-2 
(a, - Gg) x 10 


c. Figure 6.29 Force on the Pushrod F, x 10? 


1 
d. Figure 6.30 Force on the Rocker F, x 10° 
e. Figure €.31 Displacement of the Valve qa X 107! 
f. Figure 6,32 Force on the Valve F. x 10? 


i 
ge. Figure 6.33 Axial Oscillation of the Valve Stem 


i =) 
(qo < ac) Ba 8 
h. Figure 6.34 Axial Oscillations of the Pushrod 


/ -9 
(a, a a) x 10 


Figures 6,34 and 6,34 have been included for completeness and 


are not niega toe, 


138: 





The following reference points approximately correspond to 


the indicated degrees of cam rotation, Table 6.1 


Reference Degrees of 
Point Cam Rotation 


Reference Degrees of 
Poaaat Cam Rotation 





Table 6.1 Graph Reference Points and Degrees of Cam Rotation for 


Figures 6.27 - 6,32 


At point (A), the cam is at its lowest extreme and is starting 
to rise, The pushrod is peauinaatenen the cam and rocker, with 
downward motion, The valve has just undergone a severe oscillation 
and has resuned its upward motion, 

At point (B), the cam and pushrod impact and remain in contact 
for a short period, (B-C). The rising can and force of impact create 
severe oscillation of the pushrod and a bounce occurs, The pushrod 


is again thrown clear of the cam and floats (C-E). 


134 





At point (D), the valve has bounced against the valve seat, 
creating severe oscillation and more bouncing of the valve (D-E). 
Mi oushrod is floating thus is not affected by this impact, 
| At point (E), the rising pushrod now impacts against the rocker, 
and valve against the valve seat. At point (E), the valve has again 
reseated and the rocker has become almost stationary. This allows 
the rising pushrod to impact against the rocker, 

The pushrod bounces off the rocker (F-H) and again floats 
(F-H), now having downward motion and also creating downward motion 
of the valve. 

The rising cam almost immediately contacts the pushrod (G). 
The resulting motion creates positive contact between cam, pushrod, 
and rocker (H-I), and the valve bounces cease and it resumes its 
downward harmonic motion, At point (I), the speed of the cam 
throws the pushrod clear, Imnediately separation occurs at the 
rocker=— pushrod due to inertial effects, Several bounces at this 
clearance are observed (J-K-L-M), until the system reaches its 
extreme, point (M), i.e. the cam now with downward motion begins 
falling away from the pushrod, The valve spring, having reached the 
maximun extension, reverses the motion of the valve upward. The 
rocker and pushrod are now in contact (M) and both have ime 
motion. 

At point (N), the valve impacts against the seat which 
results in severe oscillations and valve bounce (P-AA). With the 


gross motion of the rocker halted, the pushrod floats clear (PAA). 


The cycle then repeats. 
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Photo 6.2 - 630 degrees Photo 6.3 - 641 degrees 





Photo 6.4 - 650 degrees Photo 6.5 - 658 degrees 
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Photo 6.6 - 664 degrees Photo 6.7 - 682 degrees 
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Photo 6.14 - 422 degrees Photo 6.15 - 432 degrees 








Photo 6.16 - 456 degrees Photo 6.17 - 468 degrees 
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Photo 6.18 - 546 degrees Photo 6.19 - 572 degrees 





Photo 6.20 - 582 degrees 
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Example 1B, Three Clearances, 11,000 RPM, No Valve Seat 


This analysis starts at 2/4 degrees of cam rotation, point 


(A). The following figures are utilized for the analysis. 


ae 


b. 


Co 
d. 
Ce 


f, 


Figume 6,35 


(1) 


(2) 


(3) 


Displacement difference between cam and pushrod 
; =2 
(a, - dg) x 10 


Displacement of lower end of pushrod, i.e. follower 


q, x 107 


Can profile, Ig x 107! 


Figure 6.36 


(1) 


(2) 


Displacement difference between pushrod and rocker 
-2 

(4, a de) x 10 

Displacement difference between cam and pushrod 


(a, - dg) x 10 


Figure 6.37 Forces on the pushrod, F, x 10° 


Figure 6.38 Forces on the rocker, F, eng 18) 


Figure 6,39 Displacement of Valve, 


; 
Z 


qe oe 10° 


Figure 6,40 Axial oscillation of the pushrod, 


(a, Gi a, ) “o 10° 


Figure 6.40 has been included for completeness and is not 


annotated, 


The following reference points approximately correspond to 


the indicated degrees of cam rotation, Table 6.2 
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Reference Degrees of Reference Degrees of 
Point Cam Rotation Point Cam Rotation 





Table 6,2 Graph Reference Points and Degrees of Cam Rotation For 
Figures 6.35~-6.39. 

The analysis is similar to example 1A, except here the motion 
of the system is not impeded by the valve seat. 

At point (A), the cam is at its lowest extreme and is starting 
to rise, The pushrod is floating between the cam and rocker, with 
downward motion, 

At point (B), the cam and pushrod impact against one another 
and remain in contact for a short period, (BC), The rising cam 
and force of impact create severe oscillations of the pushrod and 
a bounce occurs, The system remains intact, i.e. contact at the 
cam and rocker, (C-D). The valve has reached its highest extreme 
and reverses motion downwards, 

At point (D), the pushrod and cam separate. Immediately, 
contact is also lost at the rocker and the pushrod floats free, 


(D-E). 


15m 





7 eenetsae 


e 





At point (E), the rising pushrod now impacts against the rocker, 
This slows the motion of the pushrod sufficiently for the rising cam 
to again impact against the rocker, point (F). The system then 
remains in contact for a short period (F-G). 

At point (G), the speed of the cam throws the pushrod clear. 
Immediately separation also occurs at the rocker -— pushrod due to 
inertial effects. The pushrod floats free for a short time, (GH), 
and then bounces at the rocker-pushrod clearance, point (I). The 
system next reaches its extreme value, point (K), iee. the cam now 
with downward motion begins falling away from the pushrod, The 
valve spring, having reached the maximum extension, reverses the 
harmonic motion of the valve, upward. The rocker and pushrod are 
now in contact (J), and both have downward harmonic motion, 

Finally, the compression of the valve spring slows the rocker 
valve assembly until the pushrod loses contact at the rocker and 
Beain floats free, point (L). 

At point (AA), 634 degrees of cam rotation, the cycle then 


repeats. 
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Photo 6.24 - 634 degrees Photo 6.25 - 646 degrees 





Photo 6.26 - 652 degrees Photo 6.27 - 654 degrees 
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Photo 6.28 - 658 degrees Photo 6.29 - 674 degrees 





Photo 6.30 - 708 degrees Photo 6.31 - 716 degrees 
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Photo 6.32 - 364 degrees Photo 6.33 - 366 degrees 
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Fhoto 6.40 - 464 degrees Photo 6.41 - 572 degrees 
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3. Example 1C, Two Clearances, 11,000 RPM, Valve Seat 


This analysis will commence at 333 degrees of cam rotation. One 
eycle will be discussed, from point (A) to point (AA). The following 
figures are utilized for the analysis: 


a. Figure 6.41 Cam profile and displacement of pushrod, 


-l -1 
Io x LO and q, x 10 
b. Figure 6.42 
(1) Cam profile, Ig * 10 
(2) Displacement difference between cam and pushrod, 
-2 
(q) - dg) x 10 
c. Figure 6.43 Force on pushrod, F, x 10° 
d. Figure 6.44 Force on the valve, FL x 10° 
e. Figure 6.45 Displacement of valve and valve seat, 
-1 -1 
qe x LO and dio x 10 
f. Figure 6.46 Axial oscillation of the pushrod, 


(ag - a,) x 10° 


g. Figure 6.47 Axial oscillation of the valve stem, 
=3 
(q7 a) x 10 
Figures 6.46 and 6.47 have been included for completeness 
and are not annotated. 
The following reference points approximately correspond to the 


indicated degrees of cam rotation, Table 6.3. 
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Reference Degrees of Reference Degrees of 
Point Cam Rotation Point Cam Rotation 





Table 6.3 Graph Reference Points and Degrees of Cam Rotation 
For Figures 6.41 - 6.45 

The analysis is similar to that for examples 1A and 1B. For 
this example the pushrod is pinned to the rocker and thus moves as a 
part of the system rather than floating free. 

At point (A), impact against the pushrod by the cam has just 
occurred. This has created large system oscillations, first in the 
pushrod and later in the valve stem. The cam and pushrod are rising. 
At points (B) and (D) the pushrod is thrown clear of the cam due to 
its high speed of rotation. At point (E) the cam has reached its 
highest position and starts to fall. The valve has reached its lowest 
position and starts its upward travel. 

At point (F) the pushrod and rocker come into contact with very 
little force and travel together for a short distance. Separation 
again occurs at point (G) when the cam falls away from the pushrod. 
The valve in its upward travel impacts against the valve seat at 


point (H). This retards the motion of the pushrod just before contact 
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at the cam was due to occur (I). The large oscillations created 
by the valve seating cause the valve to bounce and leave contact at 
the seat, point (J). 

Meanwhile the cam has reached its lowest extreme and begins its 
upward travel. At point (K) the rising cam impacts against the 
oscillating pushrod and forces the system to return to its harmonic 


motion. The cycle then repeats, point (AA). 
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Photo 6.45 - 382 degrees Photo 6.46 - OO degrees 
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Photo 6.47 - 448 degrees Photo 6.48 - 470 degrees 
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Photo 6.49 - 482 degrees Photo 6.50 - 560 degrees 
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Photo 6.51 - 602 degrees Photo 6.52 - 604 degrees 
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Photo 6.54 - 610 degrees Photo 6.55 - 684 degrees 
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4, Example 1D, Two Clearances, 11,000 RPM, No Valve Seat 


This analysis will commence at 274 degrees of cam rotation, point 
(A). The following figures are utilized for the analysis. 


a. Figure 6.48 


al 


(1) Cam profile, ran kom 


49 
(2) Displacement of lower end of pushrod, 


q, x oa 


(3) Displacement difference between cam and pushrod, 


-2 


b. Figure 6.49, Force on pushrod, F) x 10° 
ec. Figure 6.50, Displacement of valve, 7 x on 
@. Figure 6.51, Axial oscillation of the pushrod, 
(dg - aq)) x 10~*) 
Figure 6.51 has been included for completeness and is not 
annotated. 


The following reference points approximately correspond to 


the indicated degrees of cam rotation, Table 6.4. 


Reference Degrees of Reference Degrees of 
Point Cam Rotation Point Cam Rotation 





Table 6.4 Graph Reference Points and Degrees of Cam Rotation 
For Figures 6.48 - 6.50 
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The analysis is similar to the above examples. The valve seat 
position is higher than the maximum excursion of the valve, thus 
contact does not occur between the seat and valve. The remaining 
points of interest are then only associated with clearance 1, between 
pushrod and cam. 

At point (A) the cam has reached its lowest extreme and begins 
rising. The pushrod continues its downward harmonic motion, the valve 
its harmonic motion upward. At point (B), the rising cam impacts 
against the pushrod. The motion of the ee is reversed. The 
severe oscillations and inertial effects first insure positive contact 
and then cause the pushrod and cam to separate, point (C). As the 
upward motion of the pushrod slows, the rising cam again contacts 
the pushrod, point (D). Since the motion of these two members is 
Similar, the impact force is not as severe. 

At point (E), the speed of cam rotation throws the pushrod 
off. Harmonic motion of the system then results. The cycle repeats 


at point (AA). 
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Photo 6.57 - 664 degrees Photo 6.58 - 682 degrees 





Photo 6.59 - 694 degrees Photo 6.60 - 716 degrees 
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Photo 6.61 - 400 degrees Photo 6.62 - 420 degrees 





Photo 6.63 - 430 degrees Photo 6.64 - 452 degrees 
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Photo 6.65 - 466 degrees Photo 6.67 - 586 degrees 





Photo 6.66 - 600 degrees Photo 6.68 - 626 degrees 
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5. Example two, 9000 RPM and example three, 10,000 RPM were 
choosen for comparison and completeness of this study. These examples 


are not annotated and may be found in Appendix F. 
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Vit. CONCLUSION 


The simulation for the dynamic response of an elastic link 
model with three clearances has been described. The algorithm 
developed is applied to the automobile cam actuated valve train 
assembly, with external clearances at the valve and cam, The 
internal clearance at the pushrod-left rocker arn, requires the 
solution to a second eigenvalue problem in a portion of the 
calculations. These algorithms are generally coded to aliow for 
application to similar mechanisms. The methods of calculation use 
address modification arrays to enable application of the program to 
more complicated models by changing a minimum number of variables 
and array elements. The numerous subroutines developed have, in 
general, very specific functions, and are intended to be used as 
separate calculation packages in a variety of other applications. 

While the accuracy of the mathematical model of a mechanism 
has been greatly increased over the rigid version, the complexity 
of solution has also increased, even though a number of simplifying 
assumptions were made. The major assumption made by the theory was 
that of a linear model so as to permit superposition. Other 
assumptions, made to simplify the equations but not the theory, 
involved the choice of coordinates and the use of straight, constant 
area links. 

The purpose of this work was to demonstrate a more accurate 
model of a mechanism. The extent of this work was to (1) restructure 


and reprogram an existing simulation to provide for a generalized 
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method of establishing the physical model contigurations; (2) Modity 
this simulation to allow tor more than two clearances; (3) Establish 
a method for the dynamic graphical display of the model showing the 
resulting elastic deformations of the model members at high operational 
speeds. Objectives (1) and (2) are accomplished using the IBM-360 
computer for a three clearance simulation, Dynamic dimensioning 
and program address modifications arrays contribute. to the program 
generalization. Sample problem graphical output is provided, 
Objective (3) is accomplished with the aid of the ADAGE Graphical 
Terminal, model 10 interfaced with a Scientific Data System Computer 
model 9300, Data generated by the IBM computer is tranferred to 
the SDS computer by punched cards, A program then calculates the 
model's geometry, position, and deformations. The results are then 
displayed on the AGT, These results are provided on video tape 
and may be viewed by requesting NPGS Video Tape 12NPS72000H42-VR from 
the Educational Media Department, Naval Postgraduate School. 

It is suggested that additional study be done on the calculation 
Of the pushrod bounce at the rocker. The simulation, as written, 
has difficulty in obtaining a stable system when contact is regained 
Portbhiis point, 

An extension of the AGT system would be to allow for the display 
of force and displacement curves on the AGT. Methods similar to 
those developed by this author and Lt. Bowden (Ref. 19) could be 


conveniently applied, 
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APPENDIX A 


MASS AND STIFFNESS MATRICES 


The system mass and stiffness matrices are given by equations (3.6) 


and (3.7), and are repeated below. 


9 - 
i] =) cel; tm,1, (6), : (3.6) 
i=1 
2 dl 
[K] = [Bie TE, 1 abbey (3.7) 
i=1 


The transformantion matrix [8] is given in equation (3.4) and will 
not be repeated here, The element stiffness and mass matrices are 


given by reference (1). For the axial members (pushrod, spring, and 


E, A, co 

ed 7: + 4 (A,1) 
Mj 2 -1 

[M,J, == i 7 (i=1,4,5) (A.2) 


The element stiffness and mass matrices for the bending members 


valve) they are: 





(left and right arms) are: 





i 
RE. I. 
(KJ - ii 12 -6L. 
= Te (A.3) 
a 2 
-6L, 4. 
wenn . ame 
i Al 











For equations A,1 through A.4, E, is the element elastic modulus, 
A; is the element cross-sectional area, L, is the element length, 
I. is the element moment of inertia about ote bending axis, and M, 
is the total element mass of the ‘ae member, 

Evaluating equation (3.6) for the five elements of the system 
yields the full system mass matrix, designated as layer one, The 
elements of this matrix are given by Figure A,1. Similarly, evaluating 
equation (3.7) yields the system stiffness matrix, The elements 


of this matrix are shown in Figure A.2, 


191 


























ally 

720 2 : 720 
4Lstt, : ; ~3LM, 
720 720 
156M, M ~221_M 43L.M 
» 4 5 DD 
¢ 420 * 420 720 

M 

5 

2 
~22L_M LM 3M 
720 720 720 

D 
“13L,M, -3L5M, 13M 31M ALM, 
720 720 720 720 





Figure A.1 Full System Mass Matrix 
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Full System Stiffness Matrix 
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The mass and stiffness matrices are valid for the eight co- 
ordinate system, i.e. configuration I, II, V, and VI. When the rocker 
aS in contact with the pushrod, the seven coordinate system is used 
and the transformation, mass and stiffness revert to those developed 
by Anderson, reference (2). The reduced system mass and stiffness 
matrices are shown in Figures A.3 and A.4, respectfully. These 


reduced matrices are stored in layer 3 and are used for configurations 


Mt, IV, VII, and VIII. 
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=a O @) O @) 
Oe, 22M L, ; ; -13L 
4 420 420 
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Figure A.4 Reduced System Mass Matrix 
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Figure A.4 Reduced System Stiffness Matrix 





























APPENDIX B 


NEWMARK'S 8B PARAMETER METHOD 


The one step algorithm developed by Newmark (Ref. 4), is chosen 
as the numerical integration technique used in the simulation. The 
B parameter may be changed for an analysis. This time marching 
technique solves the governing equations of motion by using the 


following relations: 


ican 


ies = (fa)? + ga) (B.1) 


and 


‘cag manera (on) ‘eu 


(B.2) 


+B (ar) (tay"** -{4") 


The net effect of the parameter B is to change the form of the 
variation of acceleration during the time interval AT. For this 
mudy, B = 8/6) assumes that the acceleration changes linearly from 


ee oe +] . 
fay? to {g¢}""-, Figure B.1. 


© 
e@ 





tt tet: 


Figure B.1 Time Step Acceleration Curve 
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In this regard equations B.1 and B.2 are written as: 


petit ees. 1 
, Sc re (a5 - a) 
Acceleration = q, (t) a a 16 (B.3) 
eyes ee 
Vel ity = 4.(t) =q.+47t _ a 4)? (B.4) 
eo AGN = ag AT 2 
+] n 
at at) 3 
‘ 7 . 1 ga eon 2 a i’ t 
Displacement = a, (t) =q,+q@,t+ 4,5 + Re ra (B.5) 
Evaluate q, (t) and q; (t) atstame t = AT, 
PP og oN «en AT , eentl AT 
, ~%TUBD*% B (B.6) 
ntl _ Ny Mam 4 en (at)* ,ogntl (AT)* (B.7) 
a PS 4! nL 
Let oA) - qd; ts qs a (B.8) 
= ye ha gam 48h (ar)* (B.9) 
37 Ty 7 Fy q; "3 
Now qn" and a can be rewritten: 
Schall n , eentl AT 
qd, =o, + 45 5 (B.10) 
ntl _ yn , ent) (aT)® eer 
il a ae 
Substitute the above results for ane and cone into the equations 
of motion 
(M] fq} +c {qa} + [kK] {q} = {F} (B.12) 


LOT 





2 
wy cay? + Cod fo + SE toy egy + Oc) cy + GEE peagyy™ 


(aie 


= {P (B.13) 


Rearranged , [to] + = (c] + or (1 | qq = ee - {cl{o}" - [Ki{y}" 


(B.14) 

rewritten, (My? cyptt = ey (B.15) 
where {P}"tt - (F} evaluated at time n+l 

[My"** = (M1 +S [cl + ar [K] (B.16) 

fae aie pe poe NT (B.17) 


This results in the accelerations fee expressed in terms of 


the known system values at the prior time step, t = 4", 


ar? .[aee"] opt 
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APPENDIX C 


FOUR CLEARANCE PROBLEM 


During the early stages of this study, limited research was 
conducted on the four clearance problem to better understand the 
necessary requirements of program generalization procedures. This 
appendix will present the unique four clearance number scheme that 
must be used in conjunction with the shifting algorithm presented 
in section 3.D. Figure ©,1 indicates the configuration number 
associated with the various combinations of clearances, for a given 
time, that the model is physically permitted to take, Figure C,2 
shows the resulting flow diagram. The shift lines relate the flow 
for a nonpositive force at the clearance indicated, Non-positive 
displacement differences are along the same shift lines, but in 


the direction opposite the force arrow. 
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APPENDIX D 


COMPUTER CONTROL CARDS 


A, IBM-360 COMPUTER CONTROL CARDS 
As previously noted subroutines were complied and placed ona 
2514 Disk Storage Device. This accomplished two goals. First, it 
eliminated the requirement to physically read into the computer a 
large number of punched cards. Second, it significantly reduced 
program execution time by greatly reducing the FORT step. This 
allowed the program to execute in a higher Ast ellacs. 
The following figures list the control cards required to accomplish 
various disk manipulations. A complete explaination will not be 
attempted here, but merely a listing of procedures that have been 
found to be effective. Section 4.4 of the Users Manual, reference 
(13), should be reviewed for more information. 
In the following example control cards, the data set name is 
90274. NILKA and is placed on public disk Mary. 
1. Initialize disk space, Figure D.1l. 
2. Load subroutine NEWMRK onto disk, Figure D.2. 
4. List disk data set members and space remaining, Figure D.3. 
4, Execute program reading only control cards, main, and data 
cards, Figure D.4. 

5. Load main program, CWG, onto disk, Figure D.5. 

6. Execute program reading only control cards and data cards, 
Figure D.6. 


7. Compress disk area when available space has been exceeded. 
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Read disk members onto scratch disk without intermediate 
Figure D.1 Initalize Disk Space 


Copy all members except NEWMRK to users disk, Figure D.9. 


Copy all members to users disk, Figure D.8. 
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Figure D.2 Load Subroutine onto Disk 
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Figure D.3 List Data Set VICC and PDS 





=250K 


N CARD) 
ITON.GO 


R /* CARD 


2314 


-SHR, UNIT= 


) 


MARY,DISP 
V 
S 


SER= 
RUP 
Ere 


? 
? 


O274.-NILKA,VOL 
B 
Y 


wis wade Ser DOK NSAT 


Wor 
wes 
OOK da 


ZN IHN <b Ul ey 
wae) «Cee Se OD Ge 


CO Or (9 <fe—- C_<t 


( ORANGE CARD ) 


In 

aon «I 
ae Oe | 
= | 


OOOIDOODOONMMDAOO IG hi Nt 


TOOWD RR Het et Re on™“. 
OF wO ee) eed ond een) sd | Cl) eee 
Ce zaA. OW a Wa B(aciat (cl eto We Oe ep ah 3G ea! 
ZO AdMe AS UO 


OLN ZSWHMALILU LU UW WL WU ue > 
Far I>: OODDOOCOrF Qa 
DQ NINANDDOIIDIDIDIONOOE WL 


s 0 J J) J 
¥~YLYOOVUUOUU UO 
FLL, PGR RE OL, FA Fc PMP AOE, Fl 


mt Od Pmt fed Pd med Pmt Ped ed pened Dom 


bk >to 
LL 7) CO 
«©e 
GO za 
Cot) | ly 
NN | 
~N™ | 


Figure D.4 Execute Program with Subroutines on Disk 
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Figure D.5 Load Main Program onto Disk 
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Figure D.6 Execute Entire Program from Disk 
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Figure D.7 Read Data Set Members onto Scratch Disk 
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Figure D.8 Copy all Data Set Members 
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Figure D.9 Copy Part of Data Set Members 





B. SDS-9300 COMPUTER CONTROL CARDS 

Transferring the SDS-9300 program, in Fortran source code, from 
punched cards to magnetic tape is a relatively easy task. foes 
accomplished, repeated compiling and execution of the program is made 
much easier. Figure D,10 indicates the control cards required to 


read the punched cards onto magnetic tape. 


-9300 COMPUTER BINARY 30CT CARD* 
DX, PAPEDT 


ome PROGRAM CARDS —--— 
=EOF 


aoe eNOUTING BARI CARDS --— 
TROUT ENE BARZ CARDS --— 
aoe OUOROULINE BARS CARDS --- 
Pee OULINE BARG CARDS --- 


=EOF 
=== SUBROUTINE GAuec Apmis) ——— 
_ETC. UNTIL ALL SUBROUTINES ARE READ. 


*Indicates card supplied by Computer Laboratory 


Figure D.10 Transferring FORTRAN Source Code From Cards to 
Magnetic Tape 


The computer laboratory provided routine TAPEDT enables the 
user to change the program coding with a minimum of cards, while 
the program remains on the tape. 

Figure D.11 shows the required control cards to compile, load, 


and execute the magnetic tape Fortran source program. 
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*Indicates Card Supplied by Computer Laboratory 


Figure D.1l1 Program Execution From FORTRAN Object on 
Magnetic Tape 


To place the program on magnetic tape as binary executable code, 
the control cards are as in Figure D.11 with an additional control 
card immediately after the 4 JOB card. Figure D.12 indicates the 


required control cards. 
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Figure D.12 Placing Program as Binary Code on Magnetic Tape 


Once the binary code is on magnetic it may be executed approximately 
15% faster than in uncompiled FORTRAN source form. The control cards 


to accomplish this are shown in Figure D.13 
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Figure D.14 Program Execution From Binary Code on Magnetic Tape 
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APPENDIX E 


TBM—360 COMPUTER PROGRAM DATA CARDS 


GNILKA, C. W. 1? DECEMBER THREE CLEARANCES NO. 1 


1,000 RPM VALVE SEAT 
1000. RPM 
7 NO CLEARANCE AT VALVE SEAT AND PUSHROD/ROCKER 
ES PRINT OUT MATRICS  MPRINT 
0 CALCOMP GRAPHS DESIRED 
0 AGT OUTPUT 
Qg NO PUNCHED CARDS; MAGNETIC TAPE 
052 DEGREES TO PRINT DEBUG OUTPUT OR STOP 
0 10 PRINTED OUTPUT DESIRED; NCOUNT 
10 PRINTED OUTPUT DESIRED FOR DEBUG PURPOSES — IPRINT 
10 | OFFLINE PRINTER GRAPH 
1 IDEG - - FORCED TO STAY IN CONFIG 7 
).000001 DELT- - TIME INCREMENT 


DISPL. DIFF. BETWEEN ROD & ROCKER ** Q(2)-Q(8) 
DISPL. DIFF. BETWEEN CAM & ROD ** Q(1)-Q(9) 
CAM PROFILE ** Q(9) 
DISPL. OF LOWER END OF PUSHROD ** Q(1) 
FORCE ON PUSHROD ** F(1) 
AXIAL OSCILLATIONS OF THE PUSHROD ** Q(8)-Q(1) 
DISPLACEMENT OF VALVE ** Q{(7) 
POSITION OF VALVE SEAT ** Q(10) 
FORCE ON ROCKER ** F(2) 
AXIAL OSCILLATIONS OF THE VALVE STEM 
FORCE ON VALVE ** F(7) 
CAM PROFILE AND DISPL. GF PUSHROD 
DISPL. OF VALVE AND VALVE SEAT 
DEGREES OF CAM ROTATION 
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APPENDIX F 


This appendix presents two additional examples for comparison 


with the 11,000 RPM examples presented in chapter 6.C. 


A, EXAMPLE 3, THREE CLEARANCES, 9,000 RPM, VALVE SHAT 
Approximately 800 degrees of cam rotation is indicated 

on the associated graphs. ‘These graphs only indicate the 

indicate the major system points of interest and are not otherwise 

annotated. The graphs, in order of sapere sia are: 


fe Figure F.1 


1 


a. Cam Profile, 10m 


49 
b. Displacement of lower end of pushrod 


-1 
qd, x 10 


c. Displacement ditference between cam and pushrod 
qd, - Ig x 107? 
2, Figure F.2, Force on Pushrod 


2 
BF, seenlO 


4. Figure F.3, Force on Rocker 


2 
Fy x 10 


4. Figure F.4, Force on Valve 


Ba x 10° 


5. Figure F.5, Displacement of Valve 


7 x 107! 


6. Figure F.6, Axial Oscillations of the Pushrod 
a =) 
(ag - a,) x 10 
7. Figure F.7, Axial Oscillation of the Valve Stem 
| -3 
- 8. 
(a, - a4) x 
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8. 


Figure F.8 
a» Displacement Difference between Cam and Pushrod 
— 
(a, - 4g) x 10 


b. Displacement Difference between Pushrod and Rocker 


=2 
(a, - 4g) x 10 


B. EXAMPLE 3, THREE CLEARANCES 10,000 RPM, VALVE SEAT 


This last example referred to in chapter 6.C. also presents 


approximately 800 degrees of cam rotation. These graphs only 


indicate the major system points of interest and are not otherwise 


annotated. The graphs, in order of presentation are: 


1. 


Figure F.9 
a. Cam Profile, dg x sofas 
b. Displacement of lower end of pushrod 


-1 
qa, x IL) 


ec. Displacement difference between cam and pushrod 


q, 7 Io mis 107? 
Figure F.10, Force on Pushrod 


2 
Fy ae ALG) 


Figure F.11, Force on Rocker 


2 
Fy pes 9, 


Figure F.12, Force on Valve 


Bo ae 10° 


Figure F.13, Displacement of Valve 


da 216 Le 
Figure F.14, Axial Oscillations of the Pushrod 


Cea) = 10 


Figure F.15, Axial Oscillation of the Valve Stem 
7 <3 
(a7 - a4) x 10 
res 





8, Figure F.16 
a. Displacement Difference between cam and pushrod 
3 
(a, we Ig) x 10 
b. Displacement Difference between Pushrod and Rocker 
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IBM-360 COMPUTER PROGRAM LISTING 


This listing is coded in FORTRAN IV, (Refs. 10-13), to run on 
an TEM-360/67, with 0S/360 release 18, The simulation requires 
234 K bytes of core storage for 850 degrees of cam rotation, 
including the CALCOMP graph buffer, The offline printer array 
(GRID) is not included in this core storage figure. Hxecution 
time is a function of the speed of rotation of the cam. For 
11,000 RPM with 850 degrees, execution time is about 4.3 minutes, 
The primary output parameters fill arrays for SYSOUT = B and 
sYsOUL = C, The B queue is for transferring data, via punched 
cards, for use in the S5DS-9300 computer. The C queue is for plotting 


graphs via the CALCOMP, 
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SDS=—9 300 COMPUTER PROGRAM LISTING 


This listing is coded in SDS FORTRAN IV, (Ref, 5-7), to run on 
a Scientific Data System Mode 9300 somputer interfaced with an ADAGE 
Graphics Display Terminal, Model 10 (AGT-10), Due to the limited 
core storage of this computer, overlaying techniques have been 


employed. 
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